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Abstract 

This paper concerns the imaging of a complex-valued anisotropic tensor 7 = ct -I- lloe from 
knowledge of several inter magnetic fields H where H satisfies the anisotropic Maxwell system 
on a bounded domain A C with prescribed boundary conditions on dX. We show that 7 
can be uniquely reconstructed with a loss of two derivatives from errors in the acquisition H. 
A minimum number of five well-chosen functionals guaranties a local reconstruction of 7 in 
dimension two. The explicit inversion procedure is presented in several numerical simulations, 
which demonstrate the influence of the choice boundary conditions on the stability of the 
reconstruction. This problem finds applications in the medical imaging modalities Current 
Density Imaging and Magnetic Resonance Electrical Impedance Tomography. 


1 Introduction 

The electrical properties of a biological tissue are characterized by 7 = u -|- te, where a and 
£ denote the conductivity and the permittivity. The properties that indicate conditions of 
tissues can provide important diagnostic information. Extensive studies have been made to 
produce medical images inside human body by applying currents to probe their electrical prop¬ 
erties. This technique is called Electrical Impedance Tomography (EIT). This imaging technique 
displays high contrast between health and non-healthy tissues. However, EIT uses boundary 
measurements of current-voltage data and suffers from very low resolution capabilities. This 
leads to an inverse boundary value problem (IBVP) called the Calderon’s problem, which is 
severely ill-posed and unstable [24, 26]. Eor IBVP in electrodynamics, we refer the reader to 
[7, 13, 20, 21, 23]. Moreover, well-known obstructions show that anisotropic admittivities cannot 
be uniquely reconstructed from boundary measurements, see [14, 26]. 
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Several recent medical imaging modalities, which are often called coupled-physic modalities 
or hybrid imaging modalities, aim to couple the high contrast modality with the high resolution 
modality. These new imaging techniques are two steps processes. In the hrst step, one uses 
boundary measurements to reconstruct internal functionals inside the tissue. In the second step, 
one reconstructs the electrical properties of the tissue from given internal data, which greatly 
improve the resolution of quantitative reconstructions. An incomplete list of such techniques 
includes ultrasound modulated optical tomography, magnetic resonance electrical impedance 
tomography, photo-acoustic tomography and transient elastography. We refer the reader to 
[1, 2, 4, 5, 15, 16, 17, 18] for more details. For other techniques in inverse problems, such as 
inverse scattering, we refer the reader to [12, 19, 25]. 

In this paper we are interested in the hybrid inverse problem of reconstructing (a, e) in 
the Maxwell’s system from the internal magnetic fields H. Internal magnetic fields can be 
obtained using a Magnetic Resonance Imaging (MRI) scanner, see [11] for details. The explicit 
reconstructions we propose require that all components of the magnetic field H be measured. 
This may be challenging in many practical settings as it requires a rotation of the domain 
being imaged or of the MRI scanner. The reconstruction of (cr, e) from knowledge of only some 
components of H, ideally only one component for the most practical experimental setup, is open 
at present. We assume that the above first step is done and we are given the data of the internal 
magnetic fields in the domain. 

In the isotropic case, a reconstruction for the conductivity was given in [22]. In [10], an ex¬ 
plicit reconstruction procedure was derived for an arbitrary anisotropic, complex-valued tensor 
7 = fj -|- te in the Maxwell’s equations in M^. This explicit reconstruction method requires that 
some matrices constructed from measurements satisfy appropriate conditions of linear indepen¬ 
dence. In the present work, we provide numerical simulations in two dimensions to demonstrate 
the reconstruction procedure for both smooth and rough coefficients. 

The rest of the paper is organized as follows. The main results and the reconstruction 
formulas are presented in Section 2. The reconstructibility hypothesis is proved in Section 3. 
The numerical implementations of the algorithm with synthetic data are shown in Section 4. 
Section 5 gives some concluding remarks. 

2 Statements of the main results 

2.1 Modeling of the problem 

Let A" be a simply connected, bounded domain of with smooth boundary. The smooth 
anisotropic electric permittivity, conductivity, and the constant isotropic magnetic permeability 
are respectively described by s{x), a{x) and /io, where s{x), a{x) are tensors and fiQ is a constant 
scalar, known, coefficient. We denote 'j = a + Luoe, where cu > 0 is the frequency of the 
electromagnetic wave. We assume that e{x) and cr{x) are uniformly bounded from below and 
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( 1 ) 


above, i.e., there exist constants > 1 such that for all ^ G 

inx 

in X. 

Let E = {E^,E'^y G and H £ C denote the electric and magnetic fields inside the domain 
X. Thus E and H solve the following time-harmonic Maxwell’s equations: 

f V X E + luiuqH = 0 , . 

\ V X H -'jE = 0, 

with the boundary condition 


uxE -.= viE"^ - V 2 E^ = f, on dX, 


(3) 


where u = (z^i, 1 ^ 2 ) is the exterior unit normal vector on the boundary dX. The standard well- 
posedness theory for Maxwells equations [8] states that given f £ H 2 (dX), the equation (2)-(3) 
admits a unique solution in the Sobolev space H^{X). In this paper, the notations V and V 
distinguish between the scalar and vector curl operators: 


X X E 


dE^ 

dxi 


dE^ 

dxo 


^ ^ ^ , OH dH., 

and XxH = {--—, -—). 

UX2 uXi 


(4) 


Although (2) can be reduced to a scalar Laplace equation for if, we treat it as a system. The 
reconstruction method holds for the full 3 dimensional case. In this paper, we assume that the 
conductivity a and the permeability e are independent of the third component in and give 
the numerical simulations in two dimension to validate the reconstruction method. 


2.2 Local reconstructibility condition 

We select 5 boundary conditions {/i}i<i <5 such that the corresponding electric and magnetic 
fields iLj}i<j <5 satisfy the Maxwell’s equations (2). Assuming that over a sub-domain 
Xq C X, the two electric fields Ei, E 2 satisfy the following positive condition, 

inf I det(£i, £^ 2)1 > Co > 0. (5) 

xGXq 

Thus the 3 additional solutions {£ 2 -i-j}i<j <3 can be decomposed as linear combinations in the 
basis {El, E 2 ), 

E2+j = )^iEi + \2^2 i 1 < J < 3, (6) 
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where the coefficients {A{, A 2 }i<j <3 can be computed by Cramer’s rule as follows: 


Ai = 


det(£? 2 +j,-£^ 2 ) _ det(V x // 2 +j, V x H 2 ) 


det{Ei,E 2 ) det(V X V X F2) ’ 
_ det{Ei, E 2 +j) _ det(V x Hi,V x i^ 2 +i) 


( 7 ) 


det(-Ei, E2 


det(V X V X H2 


Thus these coefficients can be computed from the available magnetic fields. The reconstruction 
procedures will make use of the matrices Zj defined by 




V X A{|V X A^ , where 1 < j < 3. 


( 8 ) 


These matrices are also uniquely determined from the known magnetic fields. Denoting the 
matrix := [V x i^i|V x H 2 \ and the skew-symmetric matrix J = we construct three 

matrices as follows, 

Mj = 1 < j < 3 , (9) 

where denotes the transpose of a matrix A and := {A + A^)/2. The calculations in the 
following section show that condition (5) and the linear independence of {T/j}i<j <3 in «S' 2 (C) 
are sufficient to guarantee local reconstruction of 7 . The required conditions, which allow us to 
set up our reconstruction formulas, are listed in the following hypotheses. The reconstructions 
are local in nature: the reconstruction of 7 at xq S X requires the knowledge of {Hj{x)}i<j<j 
for X only in the vicinity of Xq. 

Hypothesis 2.1. Given Maxwell’s equations (2) with smooth £ and a satisfying the uniform 
ellipticity conditions (1), there exist a set of illuminations {/i}i<i <5 such that the corresponding 
electric fields {Ei}i<i<^ satisfy the following conditions: 

1. infa,gXo I det(£i, E 2 )\ > cq > 0 holds on a sub-domain Xq C X, 

2. The matrices {T/j}i<j <3 constructed in (9) are linearly independent in 52(C) on Xq, where 
52 (C) denotes the space of 2x2 symmetric matrices. 

Remark 2.2. Note that both conditions in Hypothesis 2.1 can be expressed in terms of the 
measurements {Hj}j, and thus can be checked during the experiments. When the above constant 
Co is deemed too small, or the matrices Mj are not sufficiently independent, then additional 
measurements might be required. For the 3 dimensional case. Hypothesis 2.1 holds locally, under 
some smoothness assumptions on 7 , with 6 well-chosen boundary conditions. The proof is based 
on the Runge approximation, see [10] for details. 
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2.3 Reconstruction approaches and stability results 

The reconstruction approaches were presented in [10] for a 3 dimensional case. To make this 
paper self-contained, we briefly list the algorithm for the two-dimensional case. Denote M 2 (C) as 
the space of 2 x 2 matrices with inner product {A, B) := tr [A*B). We assume that Hypothesis 

2.1 holds over some Xq C X with 5 electric fields {Tlj}i<j< 5 . In particular, the matrices 
{Mj}i<j <3 constructed in (9) are linearly independent in 52 (C). We will see that the inner 
products of ( 7 “^)* with all Mj can be calculated from knowledge of {Rj}i<j< 5 - Then 7 can 
be explicitly reconstructed by least-square method. The reconstruction formulas can be found 
in Section 2.3.1. This algorithm leads to a unique and stable reconstruction and the stability 
estimate will be given in Section 2.3.2. 

2.3.1 Reconstruction algorithms 

We apply the curl operator to both sides of ( 6 ). Using the product rule, we get the following 
equation, 


V X Ei + Ei-V X Xi = X X E 2 +j, for j > 3. (10) 

i=l ,2 

Substituting Hi into Ei in the above equation, we obtain the following equation after rearranging 
terms. 


^VxXj- (j-^V X R,) = E - H 2 +,) (11) 

1 = 1,2 1 = 1,2 

Recalling the dehnition of Zj by (8), the above equation leads to 

7-' : = E - H 2 +j), (12) 

i=l,2 

where the matrix R = [V x RijV x H 2 ]. Note that Mj = {ZjH'^yy'^ and the RHS of the above 
equation are computable from the measurements, thus 7 can be explicitly reconstructed by (12) 
provided that {Mj}i<j <3 are of full rank in 52 (C). 

Remark 2.3. The reconstruction formulae is local. In practice, we add more measurements 
and get additional Mj such that {Mj}j is of full rank in 52 (C). The system (12) becomes 
overdetermined and 7 can be reconstructed by solving ( 12 ) using least-square method. 

2.3.2 Uniqueness and stability results 

The algorithm derived in the above section leads to a unique and stable reconstruction in the 
sense of the following theorem: 
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Theorem 2.4. Suppose that Hypotheses 2.1 hold over some Xq C X for two sets of electric fields 
{Ei}i<i <5 and {£^'}i<j< 5 , which solve the Maxwell’s equations (2) with the complex tensors 7 
and 7 ' satisfying the uniform ellipticity condition (1). Then 7 can he uniquely reconstructed in 
Xq with the following stability estimate, 

5 

II 7 - 7llvK“.°°(Xo) - ^ X/ (13) 

i=l 

for any integer s > 0 and some constant C = C{s). 

Proof. The above estimate is straightforward by noticing that two derivatives are taken in the 
reconstruction procedure for 7 . □ 


3 Fulfilling Hypothesis 2.1 

In this section, we assume that 7 is a diagonalizable constant tensor. We will take special 
CGO-like solutions of the Maxwell’s equations ( 2 ) and demonstrate that Hypothesis 2.1 can be 
fulfilled with these solutions. By definition of the curl operator, it suffices to show that 

Mj = 1 < i < 3, (14) 


are linearly independent in 52 (C), where Zj = 
following Helmholtz-type equation from (2), 


VAi|VA^2 


and H = [VHi\VH 2 ]. We derive the 


—V • 7 ^XHi + Hi = 0, for 1 < i < 5, 


(15) 


where 7 = —LLOfij'^'yJ and admits a decomposition 7 = QQ^ with Q invertible. We take special 
CGO-like solutions of the form 


Hi = 

where the Ui are vectors of unit length. Obviously, Ui defined in (16) satisfy (15) and 

0 


H = QU 


cX-Qui 

0 


oX-QU2 


(16) 


(17) 


where U = [tti|rt 2 ]- Therefore, Hypothesis 2.1.1 can be easily fulfilled by choosing independent 
unit vectors ui = ei, U 2 = 62 . Using the corresponding additional electric fields {E 2 +j}i<j< 3 , 
Cramer’s rule as in (7) yields the decompositions 

E 2 +j = X{Ei -|-A^£J 2 ) with A{ = U 2 ), A^ = ri 2 +j). 
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Then by definition of Zj, we get the following expression, 


det(u2+j, U2)(«2+j - ui), 


det{ui,U2+j){u2+j - U2)]. 


Hi '' H2 

Together with (17), straightforward calculations lead to 

ZjH'^ = H2+jQ[det{u2+j, U2){u2+j - ui),det{ui,U2+j){u2+j - U2]Q'^■ 
Using the fact that U 2 +j = {u 2 +j ■ ui)ui + {u 2 +j ■ U 2 )u 2 , the above equation leads to 


Mj — H2+jQ 


{u2+j • Ui){{u2+j ■ Ui) - 1) {u2+j ■ Ui){u2+j ' U2) 

{u2+j ■ Ul)(n2+j • U2) {u2+j ■ U2){iu2+j ' U2) - 1) 




(18) 


(19) 


( 20 ) 


where ui = ei, n 2 = 62 - Therefore, it is easy to find U 2 +j vectors of unit length such that Mj 
are linearly independent in 52 (C). 

Remark 3.1. To derive local reconstruction formulas for more general tensors (e.g. C^'°‘{X)), 
we need loeal independenee conditions of {Mj}j and we need to eontrol the local behavior of 
solutions by well-ehosen boundary eonditions. This is done by means of a Runge approximation. 
For details, we refer the reader to [3], [6] and [10]. 


4 Numerical experiments 

In this section we present some numerical simulations based on synthetic data to validate the 
reconstruction algorithms from the previous section. 

4.1 Preliminary 

We decompose 7 = cr + Love into the following form with six unknown coefficients {crj}i<j< 3 , 
{ei}i<j <3 respectively for a and e, 


7 = 


0-1 

0-2 

-\- LU ) 

' £1 

£2 ’ 

CT2 

<73 


. ^2 

^3 _ 


( 21 ) 


where each coefficient can be explicitly reconstructed by solving the overdetermined linear system 
( 12 ) using least-square method. 

In the numerical experiments below, we take the domain of reconstruction to be the square 
X = [—1,1]^ and use the notation x = (x, y). We use aN + lxN + 1 square grid with N = 80, 
the tensor product of the equi-spaced subdivision x = — 1 : h : 1 with h = 2/N. The synthetic 
data H are generated by solving the Maxwell’s equations (2) for known conductivity a and 
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electric permittivity e, using a finite difference method implemented with MatLab. We refer to 
these data as the ’’noiseless” data. To simulate noisy data, the internal magnetic fields H are 
perturbed by adding Gaussian random matrices with zero means. The standard derivations a 
are chosen to be 0.1% of the average value of \H\. 

We use the relative error to measure the quality of the reconstructions. This error is 
defined as the L^-norm of the difference between the reconstructed coefficient and the true 
coefficient, divided by the L^-norm of the true coefficient. £^., £^ with 1 < z < 3 

denote respectively the relative error in the reconstructions from clean and noisy data for <7* 
and Ei- 

Regularization procedure. We use a total variation method as the denoising procedure by 
minimizing the following functional, 

o(f) = i||f-frc||i + p||rf||Tv (22) 

where frc denotes the explicit reconstructions of the coefficients of a and e, T denotes discretized 
version of the gradient operator. We choose the /^-norm as the regularization TV norm for 
discontinuous, piecewise constant, coefficients. In this case, the minimization problem can be 
solved using the split Bregman method presented in [9]. To recover smooth coefficients, we 
minimize the following least square problem with the f^-norm for the regularization term, 

o(f) = i||f-frc||i + p||rf||i, (23) 

where the Tikhonov regularization functional admits an explicit solution f = (I + /7r*r)“^frc. 

The regularization methods are used when the data are differentiated. 

4.2 Simulation results 

Simulation 1. In the first experiment, we intend to reconstruct the smooth coefficients {uj, ej}i<j <3 
defined in (21). The coefficients are given by, 

{ (Ti = 2 + sin(7rx) sin(7ry) 

(72 = 0.5 sin(27rx) 

CTg = 1.8 + — g-15((a:+0.4)2 + (y+0.6)2) 

and 


£1 = 2 — sin(7rx) sin(7ry) 

£2 = 0.5 sin(27ry) 

gg = 1.8 _|_ e-12(a:2+y2) ^-i2{(x+0.6)^+{y-0.5)'^) 
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_ g-12((x-0.4)2 + (y+0.6)2)_ 


We performed two sets of reconstructions using clean and noisy synthetic data respectively. 
The / 2 -regularization procedure is used in this simulation. For the noisy data, the noise level is 
a = 0.1%. The results of the numerical experiment are shown in Figure 1 and Figure 2. The 
relative errors in the reconstructions are = 0.3%, = 5.1%, = 0.8%, = 33.4%, 

Eg = 0.2%, Eg = 4.9%; Eg^ = 0.1%, Eg^ = 5.8%, Eg^ = 0.5%, Eg^ = 30.0%, Eg^ = 0.1%, 
Eg^ = 4.8%. 

Simulation 2. In this experiment, we attempt to reconstruct piecewise constant coefficients. 
Reconstructions with both noiseless and noisy data are performed with /i-regularization using 
the split Bregman iteration method. The noise level is a = 0.1%. The results of the numerical 
experiment are shown in Figure 3 and 4. From the figures, we observe that the singularities of 
the coefficients create minor artifacts on the reconstructions and the error in the reconstruction 
is larger at the discontinuities than in the rest of the domain. The relative L?' errors in the 
reconstructions are Eg^ = 4.0%, Eg^ = 17.6%, Eg^ = 12.8%, Eg^ = 48.1%, Eg^ = 4.5%, Eg^ = 
16.5%; Eg^ = 0.1%, Eg = 16.3%, Eg^ = 0.5%, Eg = 35.2%, Eg^ = 0.1%, Eg^ = 16.2%. 

5 Conclusion 

We presented in this paper the reconstruction of (u, e) from knowledge of several magnetic fields 
Hj, where the measurements Hj solve the Maxwell’s equations (2) with prescribed illuminations 
/ = fj on dX. 

The reconstruction algorithms rely heavily on the linear independence of electric fields and 
the families of {Mj}j constructed in Hypothesis 2.1. These linear independence conditions can 
be checked by the available magnetic fields {Hj}j and additional measurements could be added 
if necessary. This method was used in the numerical simulations. We proved in Section 3 
that these linear independence conditions could be satisfied by constructing CGO-like solutions 
for constant tensors. In fact, these conditions can be verified numerically for a large class of 
illuminations and more general tensors. The numerical simulations illustrate that both smooth 
and rough coefficients could be well reconstructed, assuming that the interior magnetic fields 
Hj are accurate enough. However, the reconstructions are very sensitive to the additional noise 
on the functionals Hj. This fact is consistent with the stability estimate (with the loss of two 
derivatives from the measurements to the reconstructed quantities) in Theorem 2.4. 
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(i) true as 


(j) 0-3 (a = 0%) 


(k) as (a — 0.1%) 


(1) as at {y = 0} 


Figure 1: a in Simulation 1. (a)&(e)&:(i): true values of (cti, <72 , cts). (b)&:(f)&:(j): recon¬ 

structions with noiseless data. (c)&(g)&:(k): reconstructions with noisy data(Q: = 0.1%). 
(d)&:(h)&:(l): cross sections along {y = —0.5}. 
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(e) true £2 


(f) £2 {a = 0%) 


(g) £2 {a = 0.1%) (h) £2 at {y = -0.5} 



Figure 2; e in Simulation 1. (a)&:(e)&:(i): true values of (ei,£ 2 ,£ 3 )- (b)&:(f)&:(j): reconstructions 
with noiseless data. (c)&(g)&(k): reconstructions with noisy data(a = 0.1%). (d)&(h)&(l): 
cross sections along {y = —0.5}. 
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1 1.2 1.4 1.6 1.8 2 

(a) true ai 



-0.5 0 0.5 


(e) true fj 2 



1 1.2 1.4 1.6 1.8 2 


(i) true 0-3 



1 1.2 1.4 1.6 1.8 2 


(j) 0-3 (a = 0 %) 



1 1.2 1.4 1.6 1.8 2 

(c) a (a = 0 . 1 %) 



-0.5 0 0.5 


(g) (T 2 (a = 0 . 1 %) 



1 1.2 1.4 1.6 1.8 2 


(k) CT 3 (a = 0 . 1 %) 



(d) ai at {y = —0.5} 



(h) a 2 at {y = —0.5} 



( 1 ) as at {y = 0 } 


Figure 3: a in Simulation 2. (a)&(e)&:(i): true values of (cti, < 72 , cts). (b)&:(f)&:(j): recon¬ 

structions with noiseless data. (c)&(g)&:(k): reconstructions with noisy data(Q: = 0.1%). 
(d)&:(h)&:(l): cross sections along {y = —0.5}. 
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1 1.2 1.4 1.6 1.8 2 


(a) true £i 



D.5 0 0.5 


(e) true £2 



-0.5 0 0.5 


(f) £2 (a = 0 %) 



1 1.2 1.4 1.6 1.8 2 


(c) £1 {a = 0 . 1 %) 



-0.5 0 0.5 


(g) £2 {a = 0 . 1 %) 




true £2 

^-£2(0=0) 

ii 

1 

r 

1 

;j 


-1 -0.5 0 0.5 1 


(h) £2 at {y = —0.5} 



1 1.2 1.4 1.6 1.8 2 

(i) true £3 



1 1.2 1.4 1.6 1.8 2 

(j) £3 (q = 0%) 



Figure 4: e in Simulation 2. (a)&:(e)&:(i): true values of (ei,£ 2 ,£ 3 )- (b)&:(f)&:(j): reconstructions 
with noiseless data. (c)&(g)&(k): reconstructions with noisy data(a = 0.1%). (d)&(h)&(l): 
cross sections along {y = —0.5}. 
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